Supersolid Order from Disorder: Hard-Core Bosons on the Triangular Lattice 
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We study the interplay of Mott localization, geometric frustration, and superfluidity for hard-core 
bosons with nearest-neighbor repulsion on the triangular lattice. For this model at half-filling, we 
demonstrate that superfluidity survives for arbitrarily large repulsion, and that diagonal solid order 
emerges in the strongly correlated regime from an order-by-disorder mechanism. This is thus an 
unusual example of a stable supersolid phase of hard-core lattice bosons at a commensurate filling. 
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Introduction: Recent experiments on 4 He under pres- 
sure suggest Q that a supersolid phase, in which long- 
range diagonal (crystalline) order and long-range off- 
diagonal (superfluid) order coexist, may arise at temper- 
atures T < 100 mK. Such a phase had been theoretically 
envisioned long ago 0], based on the proposition that 
a nonvanishing density of zero point defects in the solid 
(intcrstitials or vacancies) could Bosc condense and form 
a superfluid "on top" of the existing crystalline back- 
ground. However, early estimates [2| of the superfluid 
fraction in clean 4 He solid were very small (~ 10 -4 ), and 
recent numerics with realistic interatomic potentials find 
no evidence for zero point defects or appreciable inter- 
particle exchanges Q in the crystal. 

Since the microscopic conditions under which clean 
He might exist in a supersolid phase are still unclear, 
it is useful to focus on simple lattice models of bosons in 
order to understand the different mechanisms by which 
a supersolid might emerge. Such "lattice supersolids" 
are superfluids which also break the lattice translation 
symmetries . Numerical studies of interacting bosons 
on the square lattice find that stable supersolid phases 
form upon doping away from a half-filled checkerboard 
or striped solid [g, 0, 0] ■ These appear to be examples 
of the "defect condensation" mechanism, where doped 
bosons (holes) act as interstitial (vacancy) sites in the 
crystal. 

In this letter we explore a different route to super- 
solidity based on the competition between Mott local- 
ization physics and geometric frustration. To illustrate 
this, we focus on a simple model of hard-core bosons 
at half-filling on the triangular lattice, interacting via a 
nearest-neighbor repulsive term, 
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For this model, we show that a supersolid phase emerges 
for V/t > 10. The supersolid phase arises from an order- 
by-disordcr effect in a strongly correlated superfluid. 
This route to supersolidity is thus complementary to the 
"defect condensation" mechanism which starts from a 
perfect crystal and considers a small defect-density in- 
duced supersolid. Our results are based on a combination 
of analytical arguments, projected boson wavefunctions, 
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FIG. 1: (color online), (a) Ground state superfluid density, 
calculated with exact diagonalization, quantum Monte Carlo, 
and variational wavefunction, with J± = 1. The inset is p s 
in the limit J z —> oo, calculated by ED (for several system 
sizes N), and VWF (in the thermodynamic limit), (b) The 
diagonal order parameter squared. Examples of the QMC 
extrapolations to N — > oo are illustrated in Fig. [5] 

quantum Monte Carlo and exact diagonalization studies. 
We thus go well beyond the earlier spin-wave work || 
on this model. It would be interesting to look for signa- 
tures of this supersolid in experiments, e.g on interacting 
cold bosonic atoms in an optical lattice [10(, where both 
a zero momentum condensate peak and Bragg peaks at 
the ordering wavevectors would be expected in the usual 
time-of-flight momentum spectroscopy. 

Alternatively, the Hamiltonian we study can be 
viewed, via the standard mapping from hardcore bosons 
to S=l/2 spins, as an XXZ model: 

H = J2 [~J± {S?S* + S?S«) + J z SfS*] , (2) 

(ij) 

with J z = V and J± = It. A spin Hamiltonian such 
as Eq. J2J may be realized in a Mott insulating phase of 
cold bosons on an optical lattice, if the bosons possess 
two internal states which can play the role of the spin 
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[lOj. Since the two-body exchange interaction for spin- 
ful bosons is naturally ferromagnetic, J± > in Eq. J5J. 
In spin language, supersolid ordering corresponds to the 
spins having their ccy-components aligned ferromagneti- 
cally (superfluid) with the z-component also ordered at a 
nonzero wavevector (solid) . In our analysis, we will work 
interchangeably in terms of bosons or spin variables. 

Spin wave theory and Landau- Ginzburg- Wilson ap- 
proach for J±/J z » 1: For large J±/J z , the spins 
align in-plane to form an XY-ferromagnetic phase. Ig- 
noring small quantum fluctuations, the ground state has 
all spins uniformly polarized along the (say) x-direction, 
giving a non-zero (S x ). This corresponds to a Bose con- 
densate of the bosons, with off-diagonal long range order 
and a nonzero superfluid stiffness. Doing a spin wave ex- 
pansion about this ordered state, we find that the gap to 
excitations at wavevectors ±ko = ±(47r/3,0) (Brillouin 
zone corners) closes at J z /J± = 2. In the boson lan- 
guage, the "roton minimum" of the superfluid hits zero 
energy. The uniform superfluid develops an instability to- 
ward 3-sublattice ordering, Sf ~ Re(V>e lk °' ri ), with the 
order parameter ip describing the critical modes. The 
Landau-Ginzburg- Wilson (LGW) expansion of the effec- 
tive action in ip and the uniform Ising magnetization M 
has the form: 
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+ ulV'l 6 +wRc(V> 6 ) + M 2 /(2 X ) - AMRe(V> 3 )] ,(3) 

where the constant x represents the boson compressibil- 
ity, and A is generally non-zero from symmetry consider- 
ations. Based as it is upon the spin- wave expansion, this 
LGW action describes the onset of non-zero ip within the 
superfluid, i.e. the ordered state is a supersolid (as op- 
posed to an insulating solid with a broken lattice symme- 
try). Exhaustive evidence supporting this interpretation 
will be provided below. 

The "massive" M field can be integrated out to obtain 
an effective action for if) alone, with the renormalizcd pa- 
rameters v = v — %A 2 /2, w — w — \X 2 /2. Two possible 
ordered states thereby occur, determined by the sign of 
w. For w < 0, ip 3 is purely real, and from Eq. a non- 
zero spontaneous uniform magnetization M is predicted. 
This ferrimagnetic state has a three-sublattice structure 
with (5?) = (m,m,-m') (M = (2m - m')/3 + 0), and 
requires phase separation in the canonical boson number 
ensemble. For w > 0, ip 3 is purely imaginary |l j |. and 
M = 0. The corresponding three-sublattice pattern has 
(Sf) = (m, —to, 0), so we refer to this as an antifcrromag- 
netic state (no phase separation is implied). Note that 
the instability to three-sublattice ordering was also found 
by Murthy et al. @ who carried out a study of the fluc- 
tuations around such possible ordered states within spin 
wave theory, and suggested a supersolid phase for large-S. 
Here we provide a more thorough and rigorous analysis, 
directly for the 5=1/2 case of most interest. Before 
proceeding, we note that Eq. © is an appropriate start- 
ing point for investigating the quantum critical behavior 



at the zero temperature superfluid-supersolid transition, 
which is expected to be of the three-dimensional XY uni- 
versality class. 

Variational wavefunction: A simple variational wave- 
function argument indicates that superfluidity survives 
in our model for arbitrarily small Jj_/J z . A good varia- 
tional approximation to the ground state for large J±/J z 
is the hard-core projection of the non-interacting Bosc- 
condensate wavefunction, 
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where AT& = N/2, and N is the number of sites. To 
account for inter-site correlations with increasing J z , we 
modify the variational wavefunction in Eq. J3J| through 
a Jastrow factor 
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where < g < 1 is a variational parameter, and {r} = 
(ri . . .rjvj. Using a Monte Carlo method, we have eval- 
uated correlation functions in ^ with the optimal g = g* 
at each value of J z /J±. 

For J z = 0, we find g* ~ 0.2, and the wavefunction 
in Eq. © has off-diagonal long range order, with (pi) ~ 
0.48. Imposing a phase gradient <f> in the wavefunction ^ 
through a twist in the boundary condition, the superfluid 
stiffness is defined as 
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where AE is the resulting change in energy. Using this 
definition, we estimate p s (J z = 0) « 0.54Ji. With in- 
creasing J z , we find g* increases and the superfluid den- 
sity decreases, as illustrated in Fig. QJi. 

Let us go to the most interesting limit of Jj_/ J z — ► 0. 
For J± = 0, the XXZ model in Eq. reduces to the clas- 
sical Ising antifcrromagnet on a triangular lattice which is 
well-known to have a macroscopic degeneracy of ground 
states « 1.381^. These are configurations with exactly 
one frustrated bond per triangle |l2j. For J z — > oo, we 
find the optimal g* — ► 1, and the wavefunction ^({r}, 0) 
is then an equal weight linear superposition of all Ising 
ground state configurations with S^ ot = 0. In this case, 
we find the off-diagonal long range order (hi) « 0.13, and 
superfluid stiffness p s « 0.09Jj_ (see Fig. HJi). This per- 
sistence of off-diagonal long range order on the triangular 
lattice for Jj_/J z — > is one of the most important points 
of this paper and is crucial for the discussion below. 

Remarkably, despite its simplicity, the variational 
wavefunction with g — > 1 can also be shown (analytically 
and numerically) to have power-law 1 / y/r correlations of 
Sf at the wavevectors ±ko, in addition to off-diagonal 
order, reminiscent of the long-range Ising correlations 
expected from the LGW theory. Since it somewhat un- 
derestimates the "solid" order, we expect that the varia- 
tional approximation slightly overestimates superfluidity, 
as also observed numerically. 
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Strong coupling analysis for Jj_/ J Z <^1: We now turn 
to directly consider the strong-coupling limit. A simple 
qualitative picture (which we will reinforce with a rigor- 
ous analysis a little later) can be obtained by mapping 
the XXZ model onto a transverse field Ising model: 



H = J, 



(ij) 



foes y^sf , 



(7) 



with h e ft — J±(Sf). We have assumed here that, as sug- 
gested by the variational wavefunction analysis, superflu- 
idity survives in the limit of small J±/ J z , This transverse 
field Ising model has been studied earlier [ill Il3l | in vari- 
ous contexts, using Ginzburg-Landau theory and Monte 
Carlo numerics. It is well established that this model 
has diagonal long-range order, in particular of antiferro- 
magnetic type, for small h e g/J z . This analogy suggests 
that the supcrsolid behavior of the XXZ model should 
be viewed, similarly to Eq. JJ}, as arising microscopically 
from an order- by-disorder mechanism |14j . 

One might worry that the above analysis was based on 
a mean- field treatment of the supcrfluid component of the 
order, which may be unreliable at small J±/J z - We now 
turn to a more exact argument, which strongly argues for 
a ground state with broken spatial symmetries. In partic- 
ular, we prove that for J±/J z — ► 0, the ground states are 
three-fold degenerate and transform non-trivially under 
lattice rotations. 

In this limit (J z = oo), the XXZ model reduces to 



Hoo = J Pc(S*S* + S?SV)P C , 



(8) 



where Pc is a projection operator onto the classical Ising 
antifcrromagnctic ground state manifold. Consider an 
L x L system with periodic boundary conditions, i.e. on 
the torus. We define the quantities E a > 0, a = 1,2,3, 
as the number of frustrated bonds along a straight closed 
path (loop) through L sites along the a th principal axis 
of the triangular lattice. These quantities are indepen- 
dent of translations of the loop. Moreover, it is straight- 
forward to show that they commute with Hqo in Eq. (JSJ 
(and indeed any local Hamiltonian in this manifold) . Be- 
cause each triangle has one frustrated bond, one also finds 
Ei + E 2 + E 3 = L. Under a six- fold rotation of the sys- 
tem, the E a cyclically permute. Hence a unique ground 
state with E\ = E2 = £3 that is also rotationally in- 
variant can be found only for L a multiple of three. For 
L 0(mod3), the ground state thus forms at least a 
three- fold degenerate rotational multiplct, as promised 
(it can be spanned by three states with "angular mo- 
mentum" I = 0, 1, 2 (mod 3)). 

More refined arguments indicate that this degeneracy 
is split only by terms of 0(e~ cL ) (with c independent 
of L) for < J±/J z <C 1, and L — > 00. Hence a non- 
trivial rotational multiplet remains the ground state with 
exponential accuracy in the thermodynamic limit. This 
behavior is inconsistent with a simple uniform supcrfluid, 
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FIG. 2: (color online). Finite-size scaling of QMC data for 
p s (a) and (m) 2 = S(ko)/L 2 (b). Here, J± = 1 and the 
inverse temperature is /3 = 10. Error bars (illustrated) are 
much smaller than the symbol sizes. Dotted lines show linear 
extrapolations for the five largest system sizes. The thick solid 
(blue) line is a fit to (m) 2 oc 1/L 2 , expected asymptotically 
when there is no diagonal LRO. 

which is expected to have only states with excitation en- 
ergies above the ground state scaling as a power-law of 
L. The most natural possibility consistent with this be- 
havior is broken spatial symmetry of the ground state. 
Other arguments can be presented against more "exotic" 
possibilities. Having concluded that the system must 
break lattice symmetries, and already argued from our 
variational wavefunction that superfluidity persists for 
J z — > 00 , we thus come to the conclusion that the ground 
state at large J z is a supersolid, as our numerical simu- 
lations below prove. 

Exact numerical results: To address the properties of 
the spin-1/2 XXZ model in Eq. @ we have carried out 
exact diagonalization (ED) studies, as well as Stochastic 
Series Expansion (SSE) quantum Monte Carlo (QMC) 
simulations in the grand canonical framework. Us- 
ing straightforward modifications to the original SSE 
directed-loop algorithm 01 1 the QMC is able to explore 
all regions of the finite-temperature phase diagram for 
small J z /J±. However for J z /J± > 10, the QMC is ob- 
served to experience a dynamical freezing due to the de- 
velopment of large energy barriers for local updates. In 
Fig. ^ we summarize our results for ED over the entire 
range of coupling strengths, and QMC in the weak to 
intermediate coupling regime (where we have been able 
to reliably equilibriate data). 

The superfiuid density is estimated by the stiffness, 
Eq. ©, easily measured using winding numbers in the 
QMC simulations 0,0]. For a given J z /J±, p s is ex- 
trapolated to L — > 00 using finite-size scaling on systems 
of size TV = L x L (Fig. . The results show that for 
small J z /J± the system is clearly in a robust supcrfluid 
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ground state with a large p s . The stiffness value de- 
creases with increasing J z /J±, however it remains finite 
to the largest values of J z /Jj_ studied with the QMC. 
The survival of p s at large J z j J± is in stark contrast to 
the behavior at J 2 /Jj_ = — 1, where the phase transition 
to a fully-polarized ferromagnetic phase is evident by the 
vanishing stiffness. Further, ED calculations (Fig.^t, in- 
set) strongly suggest that p s survives to J z /J± — > oo, 
although system sizes are too small to allow a reliable 
extrapolation to N — ► oo for this data. The small value 
of p s for J z — > oo indicates that the supersolid is very 
close in parameter space to an insulating state. 

In the above analytical arguments, the survival of p s at 
large J z / J± is accompanied by the development of simul- 
taneous Ising order (the supersolid phase). To address 
this prediction, we compute the static structure factor 

5(k) = i5^ eik ' C ' i_r<) W>- ( 9 ) 

This structure factor reveals the development of sharp 
Bragg peaks for moderate values of J z /J±. The peaks 
occur at the corners of the hexagonal Brillouin zone 
(±k = ±(47r/3, 0)), in complete agreement with the 
above analytical discussion. The presence of diagonal 
long range order (LRO) can be demonstrated rigorously 
through a finite-size scaling analysis of 5(ko) (Fig. Ht 1 ). 
Specifically, in a LRO state, the diagonal order parame- 
ter squared, (m) 2 = S(ko)/L 2 , should scale to a constant 
as L — > oo (with a leading order-correction varying as 
1/L [HI)- 111 contrast, in the absence of diagonal LRO, 
spin-spin correlations are short ranged, and (m) 2 should 
tend to zero as 1/L 2 in the asymptotic limit. Employ- 
ing this analysis on the QMC data for all J z /J± studied 
reveals the position of the zero-temperature phase transi- 
tion from the superfluid to supersolid state to be between 
J z /J± = 4 and J z /Jj_ = 5 (Fig. [J). 

Finally, we note that the grand canonical QMC simu- 
lations indicate no significant deviations from half-filling 
in the supersolid phase for system sizes up to L = 36. In 
the ED, the minimum-energy ground state is also found 
to lie within the (S z ) = subspacc. These observations 
suggest the occurrence of the antiferromagnetic rather 
than the ferrimagnetic supersolid, however QMC simu- 
lations on much larger system sizes may be necessary in 
order to distinguish conclusively between the two. 

Conclusions: We have presented a comprehensive 
study of the phase diagram of interacting hard-core 
bosons at half-filling (zero field XXZ model) on the tri- 
angular lattice. We have found that a subtle interplay 
between repulsion and quantum fluctuations produces a 
supersolid ground state in this system in the strongly cor- 
related regime. This model provides an unusual example 
of a supersolid at a commensurate density arising from 
an order-by-disorder mechanism rather than defect con- 
densation. This distinct mechanism completely avoids 
the instability of the defect-mediated square lattice su- 
persolid to phase separation ||. Here, the supersolid has 



boson density equal to or extremely close to 1/2, which 
is strongly favored by the maximal frustration-induced 
entropy of the corresponding classical Ising antiferromag- 
netic states. 
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